COMPUTING GOWDY SPACETIMES VIA SPECTRAL EVOLUTION 
IN FUTURE AND PAST DIRECTIONS 



PAULO AMORIM, CHRISTINE BERNARDI, 
AND 

^ PHILIPPE G. LeFLOCH 

O 

Abstract. We consider a system of nonlinear wave equations with constraints 
tliat arises from the Einstein equations of general relativity and describes the 
^ geometry of the so-called Gowdy symmetric spacetimes on T^. We introduce 

^ two numerical methods, which are based on pseudo-spectral approximation. 

The first approach relies on marching in the future time-like direction and 
toward the coordinate singularity t = 0. The second approach is designed 
from asymptotic formulas that are available near this singularity; it evolves 
the solutions in the past timelike direction from "final" data given at t = 0. 
This backward method relies a novel nonlinear transformation, which allows 
I I us to reduce the nonlinear source terms to simple quadratic products of the 

unknown variables. Numerical experiments are presented in various regimes, 
including cases where "spiky" structures are observed as the singularity is 
^ approached. The proposed backward strategy leads to a robust numerical 

{^JJJ method which allows us to accurately simulate the long-time behavior of a 

I I large class of Gowdy spacetimes. 

> 

1. Introduction 

We consider spacetimes {M,g) foliated by spatial slices that have the topology 
of the 3-torus and on which a two-parameter group of isometrics acts. As shown 
QQ by Gowdy [7] , a large class of such spacetimes can be described in suitably chosen 

coordinates (r, 9, a, (3) e M x in which the metric g reads 

ix! 

The unknown metric coefficients P, Q, A depend upon the variables r, 6, only, and 
are periodic in the spatial variable 6. 

In absence of matter fields, the Einstein equations for Gowdy spacetimes are 
equivalent to two second-order, nonlinear wave equations for the time-evolution of 
P,Q: 

Prr-e-'^Pgg = e'''{Ql-e-'-Ql), 
Qrr~e-^^Qge = -2{PrQr-e-^^PgQe), 
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together with first-order differential equations for the coefficient A: 

Xg^2{PrPg + e^PQrQe)- 



The functions P, Q can be determined first by solving (1.1 ), and, then, the function 



A can be obtained by integrating (1.2 1. Since, by assumption, A is a periodic 
function, w 
constraint 



function, we have J^^ Xg d9 — and, therefore, the functions P, Q must satisfy the 



/ (PrPe 
Jo 



(1.3) / iPrPe + e''^QrQe)de^O 



In fact, it is easily checked that if ( 1.3 1 holds at some positive time tq, then it holds 
for all positive times. 

The properties of Gowdy spacetimes have been extensively studied in recent 
years, both theoretically and numerically. We refer the reader to the works by 
Garfinkle [6], Isenberg [8], Isenberg and Moncrief O [10] and Rendall [13l[14], and 
Rendall and Weaver [TF, and the references cited therein. Further properties of 
Gowdy spacetimes have been established by Chae and Chrusciel jlj and Chrusciel 
and Lake |5] . The structure of Gowdy spacetimes was completely solved in a series 
of works by Ringstrom [I1]-[2T]. In particular, [TS] gives a mathematical analysis 
of the asymptotics in the expanding direction for Gowdy metrics, and [21] 
(building on the other papers) provides a complete mathematical description of 
the asymptotic behaviour for generic solutions as one approaches the singularity. 
Issues related to curvature blow-up are also well understood and generically t = 
corresponds to a curvature singularity. 

Our purpose in the present paper is restricted to presenting a new strategy 
for computing Gowdy spacetimes, and demontrating its interest for analyzing the 
behavior of spacetimes near the coordinate singularity. This strategy takes its roots 
in an approach developed by Rendall [14' jointly with Kichenessamy [TT which led 
to existence results for the above equations. Our proposed numerical approach 
will be based on spectral approximations whose applicability to problems arising in 
general relativity is now well-recognized; see the works by Ansorg [Ij, Bonazzola, 
Gourgoulhon, Grandclcment, and Novak [3J, and Lehner, Reula, and Tiglio ^12] 
and the references cited therein. 

One of the main issues of interest is to elucidate the qualitative behavior of 



solutions to the Einstein equations (1.1 1 and, especially, to understand whether 



these solutions blow-up in the past timelike direction r ^ +oo. By contrast, in the 
future direction r — > — oo the solution remains globally regular and the dynamics 
is better understood. We are here mainly interested in the region r S [tq, oo) and, 
without loss of generality, we take tq — 0. By introducing the rescaled time 

t-e-^e(0,l], 



the system (1.1) becomes 
(1.4) 



Ptt + ^-Pee = e'''{Qj-Ql), 
Qtt + Y- = -2(PtQt - PeQs)- 
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For these equations, we will consider the initial-value problem when data are pre- 
scribed on either the hypersurface t — (forward problem) , or the singularity t — 
(backward problem). 

Let us end this introduction with a brief discussion of formal expansions for the 



solutions of (1.1 1 and some heuristics arguments. 



According to [H [TOl [14] , as one approaches the coordinate singularity, the spatial 



derivative terms in (1.4 1 becomes negligible and P,Q, approach solutions of the 



ordinary differential equations 



(1.5) 



Pt 
t 

t 



-2PtQf 



These equations are referred to as the velocity term dominated (VTD) equations. 



Interestingly, the VTD equations (1.5 1 admit a large class of explicit solutions [6] 
and the asymptotic analysis of these solutions yields the following behavior (as 



(1.6) 



P = -v logt + ip + o{t), 
Q = q + t^'' (^ + o(l)). 



Here, v, ip, q and tp are prescribed coefficients depending on the parameter 9, and 
the function v is referred to as the asymptotic velocity. 



On the other hand, solutions to the full nonlinear equations (1.4 1 have been 



constructed [14 which have an asymptotic behavior described by (1.6) and are 
called asymptotically velocity term dominated (AVTD) solutions. 



It is important to observe the analysis leading to solutions of (1.4) satisfying 



(1.6) put drastic constrains for the range of the asymptotic velocity v. In fact, the 



function v must range within the interval (0, 1), except possibly at some exceptional 
points where the coefficients qg or ijj vanish. At these points, the solutions may 
develop spikes; we will discuss the properties of these spikes later in the present 
paper. The analysis is confirmed by numerous numerical experiments [6'. 

An outline of this paper follows. In Section [2] we introduce a pseudo-spectral 
approximation scheme for solving the forward initial- value problem (1.1 1. Then, in 
Section |3] we turn our attention to the backward initial- value problem for ( 1.4 1, and 
we introduce an adapted pseudo-spectral approximation for this problem. Finally, 
in Section [4] we present extensive numerical experiments; in particular, by running 
a backward evolution followed by a forward evolution, we can demonstrate the 
accuracy of the proposed method and investigate the behavior of Gowdy spacetimes 
near the coordinate singularity. 



2. Forward method based on pseudo-spectral approximation 

In spectral approximation methods, one expands the unknown function in an ap- 
propriate basis (polynomials, trigonometric polynomials, etc.) and one numerically 
solves the coupled system of equations satisfied by the components of the solution 
relative to the chosen basis. Here, we choose Fourier expansions, which are natural 
since we are dealing with periodic functions. 

Before we discretize the system (1.1) under consideration, we propose here to 
recast it in a first-order form, by using a nonlinear transformation of the dependent 
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variables: the proposed transformation allows us to reduce the nonlinear source 
terms to simple quadratic products of the unknown variables. Defining 



the Einstein equations ( 1 . 1 1 in these new variables read 
Vr-e-''We = Y^ -X^, 
Wr- e-^Ve = -W, 

2.1 

Xr~ e-^Yg ^ -(1 - V)X - WY, 

Yr ~ e-^Xg = ^VY + WX. 

Observe that the right-hand sides of the above equations contain both linear and 
quadratic terms. 

We rely on the proposed formulation (2.11 and introduce Fourier series expan- 
sions of the 27r-periodic, unknown functions V, W,X,Y, that is, we set 

(2.2) 



feez 

and analogously for the other unknowns. The Fourier coefficient V'^ is given by 
(2.3) f'=(t) = -^/ ViT,9)e-'''''d0 

while the 6'-derivative of V is obtained by differentiating each term of the above 
expansion: 

feGZ 



Then, the system (1.1 1 can be written as an infinite system of ordinary differential 
equations (ODE's) for the Fourier coefficients V'' ,W'' , X'' ,Y'^ , that is, 

V^^ - e-^ik W'' = {Y^ - X^)'' =: 

Y^^ - e-HkX^ = -{VYf + {WXf =: S^. 

The nonlinear source-terms above are (at most) quadratic in the unknowns and, 
therefore, they are easily be computed via the relation 

At this stage, we obtain the (exact!) expressions 

S'l = ^(r'rfc-' _ x^x''-^), 

iGZ 

(2.5) .q!; ^ -x'' -t-\^(v'-y''-'' -w^v'^-^) 



{AB)'' := ^A'^-'S' 



lei 

sl = ^(-T/'r'^-' -I- w^x''-^). 
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In order to arrive at a practical numerical method, two further approximations 
are necessary. On one hand, only a finite number of frequencies k can be dealt 



with in practice and, therefore, we should truncate the sum in (2.2) and by fixing a 



(sufficiently large) integer K we obtain the following approximation of the function 
V 



(2.6) 



y[^i(T,0) 



|fe|<_R' 



ike 



and we proceed similarly for the other unknowns W,X,Y. This approximation 
leads to what is referred to as the spectral method. 

On the other hand, we observe that the above approach relies on exact Fourier 
coefficients of the unknown functions, and we introduce an approximation of the 



integrals in (2.3 1 based, specifically, on a quadrature formula involving the points 



= 2Kii ■ Precisely, instead of ( |2.2[ ) , we actually define the (approximate) Fourier 



coefficients as the following averages: 



(2.7) 



1 



2K + 1 



■ikOj 



\3\<K 



and, instead of (2.6 1, we set 



(2.8) 



fee 



|fe|<_R' 



We proceed similarly for the other unknowns of the system of equations under 
consideration This approach is referred to as the pseudo-spectral method. 

Note in passing that the above method may be viewed as a collocation method. 



Namely, approximating the Fourier coefficients by (2.7 1 is equivalent to approxi 



mating the function V itself by a trigonometric polynomial which coincides with V 
on the collocation points 9j. For further details and for the corresponding theory 
of convergence, we refer the reader to the textbook by Bernardi and Maday [2 . 

Now, let us fix an integer K. Given suitable initial data for our problem, that is, 
the 27r-periodic functions 6 (V, W, X, Y){0, 6) satisfying the Einstein constraints, 
we can introduce (for all k = —K, . . . , K) 

iV,W,X,Y)^o ■=iV,W,X,Y)':{0), 



as defined in (2.7 1 from the given initial data. Then, we compute for instance 



nk 
'-'1,0 



l\rk — l 



|/|,|fe-;|<K 



and proceed similarly for the other source-terms. The range in the above summation 
is determined so that only frequencies that are below K (in modulus) are taken into 
account. 

Next, using the above expressions of the source-terms, we let successively n = 
1, 2, . . . and we evolve the Fourier coefficients by numerically solving (2.4 1 with a 
standard ODE solver. For the numerical results presented in Section |4| below, we 
have implemented both a fourth-order Runge-Kutta method and the explicit Euler 
method and finally adopted the latter since the results did not differ significantly. 
The time-discretization is defined as follows. We denote by Ar the time-increment 
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and set t„ :— nAr. Writing A instead of V, W, X, Y, and B instead of W, V, Y, X, 
respectively, we can write the Euler scheme in the form (i = 1,2,3,4) 



3. Backward method based on pseudo-spectral approximation 

3.1. Motivation for the backward strategy. We now describe the new ap- 
proach proposed in this paper. In the previous section, we described a method 
which allows us to calculate the metric coefficients P, Q numerically, by evolving 
from a non-singular coordinate time towards the singularity. We refer to this ap- 
proach as the forward method. 



We now introduce an alternative strategy and consider the problem (1.1) as 



a backward problem with data prescribed on the singularity. This allows us to 



simulate the equations (1.1) (now written in the form (1.4)) in the reverse time- 
direction. 

The behavior of solutions may then be investigated by evolving either towards 
the singularity (using the forward method) or away from it (using the backward 
method). Moreover, having two strategies allows us to evaluate the consistency of 
both methods by, first, evolving away from the coordinate singularity and, then, 
using the resulting solution as an initial data for the forward method. As it will 
turn out, this backward- forward calculation provide numerical values on the initial 
slice that are close to the initial data originally prescribed on the slice, and our 
general proposal therefore provides a reliable tool to investigate the behavior of 
Gowdy spacetimes, especially near the singularity. 

It is important to stress that performing this backward-forward evolution does 
not amount to simply computing the solution from t = to t = 1 and, then, 
"undoing" the computation in the opposite direction. Indeed, the two systems of 
equations involved in the numerical simulation are highly nonlinear and, further- 
more, they evolve the solution in opposite time directions. 

Another validation of the proposed backward strategy is given by the fact that, 
as we will see, the numerical results of our forward simulations are virtually indis- 
tinguishable from the ones found in the literature. See Section |4] below for further 
details. 

In order to move away from the singularity, we start by considering the Einstein 



equations in the form (1.4). However, if we try to convert this into a first-order 



system in the same way as in the previous section, we obtain 

Vt-We^^iX^ -Y^), Wt-Ve = -W, 

Xt-Ye=^-{{l-V)X + WY), Yt-Xe=^-{VY -WX). 

This system is not amenable to stable numerical simulations as it stands, due to 
the singularities 1/t arising in the right-hand sides. Indeed, if we were applying 
the Euler method and evolve the solutions in time, the scheme would not depend 
on the time-increment. Moreover, as can be seen from the asymptotic expansions, 
three of these four quantities vanish (in general) on the coordinate singularity. 
Therefore, it is not possible to use the above formulation to numerically evolve 
from the coordinate singularity. 
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3.2. Rescaled initial-value problem. Hence a new approach, presented now, is 
needed to successfully evolve from the singularity. Recall that the behavior of the 
metric coefficients is determined by the main coefficients w, ^/j, g of the asymptotic 
expansions of P,Q. This motivate us to introduce the rescaled unknowns 



(3-1) ,l-2„ 



V ^ -tPt, ^^P + vlogt, 



which, according to the asymptotic expansions (1.6 1, are expected to have a limit 
on the singularity and to converge toward the given data w, (/s, t/i, g as t — s- 0: 

^it,9)^2v{e)m, x{t,9)^q'{9). 

Expressing the Einstein equations as a system of evolution equations in terms of 
the proposed unknowns yields the system 

Vt = -me + tv" log t - e^* (^2-1^1/2 _ t^-^-x^) , 
We then consider the initial- value problem with data prescribed on t = and, more 



precisely, given data v, V', 9, we solve (3.2 1 for t > 0, with initial data 



V{Q,9) = v{9), $(0,0) = ^(0), 

^(0, 9) = 2v{9m9). X(0, 9) = q'{9). 

3.3. Heuristics revisited. The condition v S (0, 1) arises formally from the for- 



mulation (3.2 1 of the Einstein equations. Indeed, we expect V to have a finite limit. 



V, on the slice i = 0. However, the first equation in (3.2) is of the form Vt ^ t" , 
with V — min(2w — 1, 1 — 2v), unless ^' or x vanish on the singularity. This suggests 
(at least formally) that V ^ i""*"^ -I- Vq, and according to this, V has a finite limit 
on i = if and only if u e (0, 1). 

Thus, if V has a finite limit on the coordinate singularity, then this limit must 
lie in the interval (0, 1), except possibly on a set of exceptional values where either 
^' or X vanish. 

Similarly, in view of the third and fourth equations we see that if (for instance) 
?; > 1 at some point, then xe must vanish at that point -in order to ensure that 
^' has a finite limit on the singularity. Similarly, must vanish where f < 0, to 
ensure that x has a finite limit for t — 0. 



We thus see that, at least formally, the unknowns in (3.2 1 have finite limits on 
the coordinate singularity if either ?; € (0, 1) or 

V > I and X — Xe 0, 
t; < and ^^ig = 0. 



We will see that our numerical results fully agree with this heuristical analysis. 
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3.4. Pseudo-spectral scheme. We approximate the Einstein system for the Gow- 



dy spacetimes, ( 3.2 ) , ( 3.3 ) by using a Fourier pseudo-spectral approximation scheme. 



We follow the main lines of the schemes developed in Section [2] however some addi- 
tional care is needed for the discretization of the (singular) equations (3.2). Inter- 



estingly enough, the heuristics discussed previously provides the necessary insight. 

Denote by At the time-increment and set i„+i/2 := {n + 1/2) At. Consider, for 
instance, the first equation in (3.2 1. Applying the Fourier transform and using the 



explicit Euler scheme for the time derivative, we arrive at the following equations 
for all k = -K, ...,K 



--v::-{n+i/2){ikAtr^^, 

-(e^*"(( 



l/2){ikAtyv''logt 



-1/2 



l/2)2"-i(At)2^«'2 - (n+ l/2)i-2''(At)2-2«..2 



The key observation here is that all powers of At are positive, at least as long as 
the velocity is chosen to satisfy v £ (0, 1). Hence, the above method makes sense 
under precisely the same conditions we derived by formal analysis. 



4. Numerical experiments 



4.1. Asymptotic velocity inside the interval (0, 1). In this section, we provide 
extensive numerical tests performed with the methods proposed in this paper, and 
we implement the backward- forward strategy sketched earlier in Section [3] 

We begin with the case that the asymptotic velocity lies in the interval (0,1). 
We run the backward method described in Section [3] with initial data 



(4.1) 



V{0) =0.2 sin 361 -I- 0.1 cos 56* -I- 0.5, 

$(0) ^ 0.35 sin 46* + 0.1 cos 16* + 0.2 sin 76* + 0.3 cos 56* + 0.5 sin 36* + C, 
^-(0) = 0.25 sin 46* + 0.5 cos 39 + 0.4 sin 29 + 0.2 cos 66* + 0.15 sin 76, 
x(O) = -sin20, 



where the constant C is chosen so that the constraint (1.3 1 holds. In Figure [T] we 



display the solution computed by the backward code at the time t — 1. 

Next, as outlined in the previous section, we use the above result as an initial 
data at T = for the forward method. Finally, at the time r = — log(A</2), 
we can compare the numerical asymptotic velocity V with the first iteration of 
the backward code. At r = — log(At/2) + 1, we can also compare the numerical 
asymptotic velocity V with the initial data originally prescribed in the backward 
code. We proceed similarly for the other unknowns. These last two results have 
been found to be completely similar, and so we plot the corresponding error at the 
time T = -log(At/2) + 1, only. See Figure [2] 

In Figures [3] and |4] we plot the maximum value of this error as a function of the 
discretization parameter K (recall that the total number of points is 2K + 1) for 
various choice of time steps. 

In Figure [5] we show the evolution of the constraint (1.3 1 during the backward 
and the forward evolutions. 



4.2. Asymptotic velocity outside the interval (0, 1). 
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Figure 1. Solution computed by the backward code at t = 1 



(r = 0), with initial data given by (4.1 1 



4.2.1. Agreement with previous numerical results. We now use our numerical meth- 
ods to investigate the exceptional situation in which the asymptotic velocity v takes 
some values outside the interval (0, 1). As observed in the heuristic discussion above, 
this may only occur at those points where either or x vanish on the singularity. 
Furthermore, to ensure that \E' and x approach to finite and regular limits, we must 
assume that X9 f^^d also vanish at those points. 

Figure |6] displays the metric coefficients P, Q computed by the forward code for 
the (large) initial data 



(4.2) 



T/(0) = 5cos6l, X(0) = -sin( 



r(0) = W{0) ^ 0, 



which are known to produce spikes in the velocity (and, hence, also in P and Q). 
Our results are in complete agreement with the previous work [B]. 



4.2.2. The case that x vanishes on an interval. We can validate the heuristic anal- 
ysis by relying on the proposed backward-forward strategy. We choose the initial 
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Figure 2. Error between the first iteration of the backward 
method and the numerical values after backward and forward evo- 
lutions, T = -log(At/2). 



data 



V{0) = 0.2 sin 6* + 0.05 sin 56* 
$(0) = cos 61, 

^'(0) = sin(6l - tt/2) + 0.76, 

'o, 



0.95, 



x(o) = 



(?e (0,7r), 
0.2sign(sin26l)sin(26')2, 6*6 (7r,27r), 

where the velocity v — V{0) is larger than 1 on some interval; see Figure [t] Ac- 
cordingly, we set X = on that interval, but we note that x does not vanish on the 
whole interval [0,27r]. The difference between these initial data and the result of 
the backward- forward simulation is plotted in Figure |8] 

It is instructive to consider also a case where the backward-forward approach 
does not accurately capture the behavior of the solutions. Namely, by modifying 
the initial velocity V{0) above, so that it is now larger than 1 in some interval where 
X does not vanish. According to our heuristics, it is not possible to obtain these 
data on the singularity by evolving towards it. Nevertheless, the backward method 
apparently does converge and generates some solution on the hypersurface t = I. 
By taking this numerical results as the initial data for the forward method and 
then evolving it toward the singularity, we have found that the computed solution 
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Figure 3. Error between the backward- forward solution and the 
first iteration of the backward code at r = — log(Ai/2), as a func- 
tion of the number of points K (for different time steps). 



ends up being very far from the initial data originally prescribed. This suggests 
that the solution computed by the backward code was physically meaningless. Our 
approach allows us to detect such a behavior. 

4.2.3. The case where vanish at a single point. The initial value v for the 

backward code is now chosen to lie outside the interval (0, 1) at one point, only. The 
purpose here is to show that our method can accurately simulate these exceptional 
situations. We perform the backward-forward simulation and compare the result 
to the prescribed velocity. 
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Figure 4. Error between the backward- forward solution and the 
initial data of the backward code at t = — log(Ai/2) + 1, as a 
function of the number of points K (for different time steps). 





Figure 5. Evolution of the constraint (1.3 1 during the backward 
evolution (left) and the forward evolution (right). 
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Figure 6. P (left) and Q (right) at t = 10, as computed by the 
forward code for initial data producing spikes. 
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Figure 7. Initial data for the backward method with v > 1 and 
X = on an interval. 



where A = 20, and the constant C is chosen so that the constraint (1.3 1 holds. 
At the point 9 = n where v equals 1, both x ^iid xg vanish, as required by the 
heuristics. The error between the computed quantities after the backward-forward 
simulation and the initial data for the backward method is shown in Figure [9] 

This initial data chosen for the velocity v is intended to qualitatively resemble 
the spiky features observed in Gowdy spacetimes, which are discussed below in 
more detail. 

It is interesting to also consider various values of the parameter A, and to look for 
a possible limit beyond which our simulations are no longer accurate. As A grows, 
the asymptotic velocity V{0; A) approaches the constant 1/2, except at 61 = tt, where 
it takes the value 1. This is thought to be the asymptotic behavior of the spikes on 
the velocity. 

It is natural that our method breaks down beyond a certain value of A, since 
the initial velocity is then approaching a singular function. Still, we have found 
that the approximation by the backward-forward method still works up until about 
A — 900 (with K — 500 and time step 0.0005). However, this accuracy is improved 
by taking larger values of K and smaller values of the time increment. 
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Figure 8. Error between Y^^^^^x given on the singularity and 
the numerical values after backward and forward evolutions, r — 
11.5, with V > \ and % = on an interval. 



The solution of the backward simulation sX t — \ for the initial data (4.3 1 and 
A = 1000 is shown in Figure [TO) The sharp gradient features observed in the 
solution at the time t = \ are simply due to the fact that the data itself was 
chosen to be have a sharp gradient on the singularity. In Figure we show the 
corresponding spiky upward-pointing initial data v (left) and the error between the 
numerical velocity after a backward-forward simulation and the initial data for the 
backward method (right). In this case, the error in the variables ^',x is rather 
large, showing that, in these more challenging cases, some information is lost in 
the backward- forward evolution. Nevertheless, this discrepancy between the initial 
data for and the values computed after the backward- forward evolution is 
localized on small regions of the interval [0,27r]. The functions ^'(0) and x(0) are 
otherwise well approximated. 

Since our main goal, at this stage, is not to provide a competitive numerical 
scheme for Gowdy spacetimes, but rather to demonstrate the validity of the pro- 
posed strategy, we argue that this discrepancy does not invalidate our strategy. It 
merely shows that our method breaks down when dealing with steep, singular data 
-which is not surprising from the standpoint of spectral approximation theory. 

The situation here is quite different from the one described in the previous para- 
graph, where we expected a large difference between, say, the data ^(0) and the 
result of the backward-forward evolution. In that case, we observed that the two 
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resulting functions were completely unrelated. In contrast, in the present situation 
our approach preserves the overall shape of the functions, at least, and is also a 
rather good approximation in most regions. 

Finally, note that the asymptotic velocity v (which is well approximated) is the 
most important quantity from the physical standpoint, since (see v is in fact 
the only coefficient which determines the nature of the coordinate singularity: the 
value of V alone determines whether the curvature blows up as one approaches the 
coordinate singularity. Therefore, even in these cases where our method does not 
provide an accurate approximation of certain coefficients, the physically relevant 
coefficient v is always well approximated, at least in the variety of cases we tested. 

The results obtained with initial data satisfying u < at some point are similar 
to the ones for u > 1 but, now, the data ^'(0) must be such that both 'i', "ifg vanish 
at those points. The error between the initial data and the computed backward- 



forward quantities is plotted in Figure |12[ for initial data similar to (4.3 1 and with 
A = 20. 
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Figure 9. Error between \E',x given on the singularity, and 
numerical values after backward and forward evolutions, r 11.5, 
for V = 1 a,t one point. 



4.2.4. Relation to spikes on the velocity. Since the points where ^' or x vanish 
are, in general, isolated points, the asymptotic velocity v generated by the forward 
code (and, in turn, the metric variable P) may exhibit spikes at these exceptional 
points. See, for instance, Figure [6] This happens because the solution V of the 
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Figure 10. Solution computed by the backward code at t = 1 



(r = 0) with initial data given by (4.3 1. (This solution is taken as 



the initial value for the forward code.) 
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Figure 11. Spiky initial asymptotic velocity, i = (left), and the 
error between the computed asymptotic velocity after backward 
and forward evolutions and the initial data, r — 10. 
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Figure 12. Error between <i>, ^, x given on the singularity, and 
numerical values after backward and forward evolutions, r = 11.5, 
for w = at one point. 



forward problem (2.11 (which tends to w as t — > 0) is forced to enter the interval 
(0, 1) everywhere except at these isolated points, where it may remain outside the 
interval. This is the origin of the spiky features observed in Gowdy spacetimes. 
These spikes have been extensively studied both numerically and theoretically. 

The main purpose of the previous simulations was to demonstrate that the 
forward-backward strategy allows one to deal with situations involving spikes, while 
still maintaining an acceptable level of accuracy. We emphasize again that the fact 
that accuracy breaks down when the spiky initial formation is sufficiently steep is 
completely natural, since our method needs sufficiently smooth data to work well, 
and that given sufficient spatial resolution and computing power, we should be able 
to accurately simulate any given spiky feature by the forward-backward strategy. 

4.2.5. Limitations and extensions of the method. As a natural extension of our 



approach we could choose some initial data for the forward method, such as (4.2 1, 
then evolve it to a large time r, and finally apply the backward method to the 
result. As estabUshed in the present paper, this strategy is successful for initial 
data originating from solutions without spikes. However, our formulation does not 
allow us to consider more general data. To see this, observe that, in Figure |6] the 
metric coefficient Q appears to develop a discontinuity near ^ = 1.5. This means 
that Qg will blow up on the coordinate singularity. On the other hand, in view 
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of the definition of the unknowns for the backward method (see (3.1 1), this means 
that X does not take a finite value on the singularity. Thus, our backward method 
does not allow us to simulate such solutions. 

In principle, it should be possible to use the forward-backward strategy and 
design a numerical method covering such cases. We could define new unknowns (in 
the spirit of ( |3.1[ ) but) taking into account higher-order terms in the corresponding 
expansions. For instance, we could introduce 

which, according to the asymptotic expansions ( |1.6[ ), should now converge to a 
finite limit, even in cases such as the one in Figure [6) 



5. Concluding remarks 

The numerical simulation of (coordinate) singularities in general relativity is a 
particularly challenging problem, from both the mathematical and the computa- 
tional standpoints. In the present paper, we have initiated a new strategy for the 
numerical investigation of Gowdy spacetimes: 

• We have recovered classical results using a Fourier pseudo-spectral scheme 
which produces accurate results using rather modest computing power. In- 
deed, these simulations were performed on a MacBook Pro at 2.4GHz, using 
the free software Scilab. The running time was of the order of 2-3 hours 
for the longest computations. 

• We have designed a new backward scheme, also based on pseudo-spectral 
approximation, which is computationally undemanding and allows to sim- 
ulate the asymptotic behavior by considering an initial-value problem with 
"final data" on the coordinate singularity. 

• This backward approach was validated by first evolving away from the 
singularity (using the backward code) and then evolving back toward the 
singularity (using the forward code). 

After a backward-forward evolution, we have obtained accurate results, and we 
firmly established the backward approach as a reliable approach to study Gowdy 
spacetimes. Interestingly enough, the proposed backward-forward strategy fails to 
work in the situations where it is not assumed to work, that is, when we set as 
"final data" on the singularity some data which are not allowed by the heuristics 
described in Section [331 

The Einstein equations for Gowdy spacetimes are known to produce solutions 
with non-smooth behavior on the coordinate singularity (cf. the spikes discussed 
in the previous section) . Through extensive numerical experiments we have estab- 
lished that the backward-forward simulation is valid even with initial (backward) 
data which are "almost spiky" , while still retaining a good quality of the results. 
Our method allows us to cover a large class of spikes, and, moreover, suggestions 
have been made on how to extend the proposed scheme and possibly deal with even 
more singular regimes. 

In short, the proposed backward strategy leads to a robust numerical method 
which allows us to accurately simulate the behavior of a large class of Gowdy 
spacetimes near the singularity. 
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